package org.hit.burkun.swalk.simulation;

import java.util.ArrayList;

import cern.colt.matrix.tdouble.DoubleMatrix1D;
import cern.colt.matrix.tdouble.impl.SparseCCDoubleMatrix2D;

/**
 * 
 * Graph implementation
 * 
 */
public class LapcNetwork extends RandomWalkGraph {
	// useful
	/** Sum of the each row of the adjacency matrix */
	public double[] rowSums;
	public SparseCCDoubleMatrix2D[] dLaplac;
	double[] sumR = new double[dim];

	/**
	 * Constructor
	 * 
	 * @param n
	 *            : number of nodes
	 * @param f
	 *            : number of features
	 * @param s
	 *            : index of starting node
	 * @param list
	 *            : the this given as an array of features
	 * @param D
	 *            : the linked set
	 * @param L
	 *            : the no-link set
	 */
	public LapcNetwork(int n, int f, int s, ArrayList<FeatureField> list,
			ArrayList<Integer> D, ArrayList<Integer> L) {
		super(n, s, f, list, D, L);
		this.rowSums = new double[this.dim]; // filled in the
												// buildTransitionMatrix method
		this.dLaplac = new SparseCCDoubleMatrix2D[f];
		for (int j = 0; j < f; j++) {
			// 计算laplc
			dLaplac[j] = new SparseCCDoubleMatrix2D(dim, dim);
		}
	}
	
	public LapcNetwork(RandomWalkGraph rg) {
		super(rg.dim, rg.s, rg.f, rg.list, rg.D, rg.L);
		this.rowSums = new double[this.dim]; // filled in the
												// buildTransitionMatrix method
		this.D = new ArrayList<Integer>();
		this.L = new ArrayList<Integer>();
		this.dLaplac = new SparseCCDoubleMatrix2D[f];
	}

	/**
	 * Build the Adjacency matrix of the graph
	 * 
	 * @param param
	 *            : the parameters used for building the adjacency matrix
	 * @return
	 */
	@Override
	public void buildAdjacencyMatrix(DoubleMatrix1D param) {
		for (int j = 0; j < f; j++) {
			// 计算laplc
			dLaplac[j] = new SparseCCDoubleMatrix2D(dim, dim);
		}
		A = new SparseCCDoubleMatrix2D(dim, dim);
		double temp;
		int r, c;
		sumR = new double[dim];
		for (int i = 0; i < list.size(); i++) {
			r = list.get(i).row;
			c = list.get(i).column;
			double value = param.zDotProduct(list.get(i).features);
			temp = weightingFunction(value);
			A.set(r, c, temp);
			sumR[r] += temp;
			if (r != c){
				A.set(c, r, temp);
				sumR[c] += temp;
			}
			//必须在A之后
			for (int j = 0; j < f; j++) {
				// 计算laplc
				dLaplac[j].set(r, c, weightingFunctionDerivative(i, r, c, j));
				if(c!=r){
					dLaplac[j].set(c, r, weightingFunctionDerivative(i, c, r, j));
				}
			}
		}
		for (int i = 0; i < list.size(); i++) {
			r = list.get(i).row;
			c = list.get(i).column;
			double norValue = A.get(r, c) / Math.sqrt(sumR[r] * sumR[c]);
			A.set(r, c, norValue);
			if (r != c)
				A.set(c, r, norValue);
		}

	}

	/**
	 * Builds the transpose of the transition matrix for given adjacency matrix
	 * and damping factor
	 * 
	 * @param alpha
	 *            : damping factor
	 * @return SparseCCDoubleMatrix2D
	 */
	@Override
	public SparseCCDoubleMatrix2D buildTransitionTranspose(double alpha) {

		SparseCCDoubleMatrix2D Q = new SparseCCDoubleMatrix2D(this.A.rows(),
				this.A.columns());

		// row sums
		int r, c;
		for (int i = 0; i < this.dim; rowSums[i++] = 0)
			;
		for (int i = 0; i < this.list.size(); i++) {
			r = this.list.get(i).row;
			c = this.list.get(i).column;
			rowSums[r] += this.A.get(r, c);
			if (r != c)
				rowSums[c] += this.A.get(c, r);
		}

		// (1-alpha) * A[i][j] / sumElements(A[i])) + 1(j == s) * alpha
		// build the transpose of Q
		double value;
		for (int i = 0; i < this.list.size(); i++) {
			r = this.list.get(i).row;
			c = this.list.get(i).column;
			value = this.A.get(r, c);
			value *= (1 - alpha);
			// value /= rowSums[r];
			// Q.set(c, r, value);
			Q.set(c, r, value / rowSums[r]);

			if (r == c)
				continue;

			// value = this.A.get(c, r);
			// value *= (1 - alpha);
			// value /= rowSums[c];
			Q.set(r, c, value / rowSums[c]);

		}

		for (int i = 0; i < Q.rows(); i++) {
			value = Q.get(this.s, i);
			value += alpha;
			Q.set(this.s, i, value);
		}
		return Q;
	}

	/**
	 * Returns matrix of partial derivatives of the transition matrix with
	 * respect to the featureIndex-th parameter for the given graph
	 * 
	 * @param featureIndex
	 *            : the index of the parameter with respect to which the
	 *            derivative is being calculated
	 * @param alpha
	 *            : the damping factor
	 * @return SparseCCDoubleMatrix2D
	 */
	@Override
	public SparseCCDoubleMatrix2D transitionDerivativeTranspose(
			int featureIndex, double alpha) {

		SparseCCDoubleMatrix2D dQt = new SparseCCDoubleMatrix2D(this.dim,
				this.dim);

		// derivative row sums
		double[] cache = new double[list.size()];
		int r, c;
		double[] dRowSums = new double[this.dim];
		for (int i = 0; i < this.list.size(); i++) {
			r = this.list.get(i).row;
			c = this.list.get(i).column;
			double weid = getLapDerivative(i, r, c, featureIndex);
			cache[i] = weid;
			dRowSums[r] += weid;
			if (r != c)
				dRowSums[c] += weid;
		}

		double value;
		for (int i = 0; i < this.list.size(); i++) {
			r = this.list.get(i).row;
			c = this.list.get(i).column;
			double weid = cache[i];
			value = (weid * rowSums[r]) - (this.A.get(r, c) * dRowSums[r]);
			value *= (1 - alpha);
			value /= Math.pow(rowSums[r], 2);
			// dQ.set(r, c, value);
			dQt.set(c, r, value);

			if (c == r)
				continue;

			value = (cache[i] * rowSums[c]) - (this.A.get(c, r) * dRowSums[c]);
			value *= (1 - alpha);
			value /= Math.pow(rowSums[c], 2);
			// dQ.set(c, r, value);
			dQt.set(r, c, value);
		}
		return dQt;
	}

	@Override
	public double weightingFunction(double x) {
		// return Math.exp(x);
		if (x > 30)
			return 1.0;
		if (x < -30)
			return 0.0;
		return 1.0 / (1.0 + Math.exp(-x));
	}

	/**
	 * Calculate partial derivative of the weight function (exponential funcion
	 * considered) parameterized by w, with respect to the index-th parameter
	 * for the given graph
	 * 
	 * @param nodeIndex
	 *            : the index of the node in the graph
	 * @param row
	 *            : the row index of the adjacency matrix
	 * @param column
	 *            : the column index of the adjacency matrix
	 * @param featureIndex
	 *            : the index of the parameter with respect to which the
	 *            derivative is being calculated
	 * @return double
	 */
	@Override
	public double weightingFunctionDerivative(int nodeIndex, int row,
			int column, int featureIndex) {
		double x = A.get(row, column);
		double featureValue = this.list.get(nodeIndex).features
				.get(featureIndex);
		if (x > -100 && x < 100) {
			double ex = Math.exp(x);
			return ex / ((1.0 + ex) * (1.0 + ex)) * (featureValue);
		} else {
			return 0.0;
		}
	}

	public double weightingFunctionDerivativeLate(double x, double featureValue) {
		if (x > -100 && x < 100) {
			double ex = Math.exp(x);
			return ex / ((1.0 + ex) * (1.0 + ex)) * (featureValue);
		} else {
			return 0.0;
		}
	}

	/**
	 * Returns true if a link from 'from' node to 'to' node in the graph,
	 * otherwise returns false
	 * 
	 * @param from
	 *            : link start node
	 * @param to
	 *            : link end node
	 * @return boolean
	 */
	@Override
	public boolean hasLink(int from, int to) {
		return this.A.get(from, to) > 0 || this.A.get(from, to) < 0;
	}

	private double getLapDerivative(int nodeIndex, int row, int column,
			int featureIndex) {
		double fenMu = sumR[row] * sumR[column];
		double weiFunValue = A.get(row, column) * Math.sqrt(fenMu);
		double sqrtSum = Math.sqrt(fenMu);
		double featureValue = this.list.get(nodeIndex).features
				.get(featureIndex);
		double dArc = weightingFunctionDerivativeLate(weiFunValue, featureValue);
		// 对（Sum(Aju)*Sum(Avi)求导
		ArrayList<Integer> ajR = new ArrayList<>();
		for (int i = 0; i < dim; i++) {
			if (i != row) {
				if (hasLink(row, i)) {
					ajR.add(i);
				}
			}
		}

		ArrayList<Integer> ajC = new ArrayList<>();
		for (int i = 0; i < dim; i++) {
			if (i != column) {
				if (hasLink(i, column)) {
					ajC.add(i);
				}
			}
		}
		double dSqrtSum = 0;
		for (int i = 0; i < ajR.size(); i++) {
			int rowIdx = ajR.get(i);
			double ridxDArc = dLaplac[featureIndex].get(row, rowIdx);
			double ridxArc = A.get(row,rowIdx) * Math.sqrt(sumR[row] * sumR[rowIdx]);
			for (int j = 0; j < ajC.size(); j++) {
				int colIdx= ajC.get(j);
				double cidxDArc = dLaplac[featureIndex].get(colIdx, column);
				double cidxArc = A.get(colIdx, column) * Math.sqrt(sumR[colIdx] * sumR[column]);
				dSqrtSum += ridxArc * cidxDArc + ridxDArc * cidxArc;
			}
		}
		dSqrtSum = dSqrtSum / (2 * (sqrtSum));
		return (dArc * sqrtSum - dSqrtSum * weiFunValue) / fenMu;

	}
	
	// TODO: DEBUG
		@Override
		public void printMatrix(SparseCCDoubleMatrix2D mat) {
			for (int i = 0; i < mat.rows(); i++) {
				for (int j = 0; j < mat.columns(); j++)
					System.out.printf("%5.4f ", mat.get(i, j));
				System.out.println();
			}
			System.out.println();
		}

		// TODO: DEBUG
		public void isColumnStochastic(SparseCCDoubleMatrix2D mat) {
			System.out.println();
			for (int i = 0; i < mat.columns(); i++)
				System.out.printf("%.2f\n", mat.viewColumn(i).zSum());
		}
}
